
R version 4.5.0 (2025-04-11) -- "How About a Twenty-Six"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

- Project '~/repos/wong_bjps_2025' loaded. [renv 1.1.5]
> ## This file assesses the matches with whites using all possible DAs in canada.
> options(digits = 4, width = 132)
> 
> library(here)
> library(tidyverse)
> library(RItools)
> 
> source(here::here("Design", "nonbimatchingfunctions.R"))
> 
> load(here::here("Design", "matches_anyDA_new.rda"), verbose = TRUE)
Loading objects:
  matches_anyDA_new
> load(here::here("Data", "wrkdatOwnMap_anyDA_new.rda"), verbose = TRUE)
Loading objects:
  wrkdatOwnMap_new
> 
> wdat0 <- wrkdatOwnMap_new %>%
+   filter(!is.na(vm_change) & !is.na(social.capital01) & !is.na(community.resp01)) %>%
+   droplevels()
> 
> wdat0 <- left_join(wdat0, matches_anyDA_new)
> table(is.na(wdat0$pair))

FALSE  TRUE 
 3772  2614 
> 
> ## Easier to think about "higher perceiver" versus "lower perceiver" when it comes to asking whether the matching worked.
> ## Further simplify the working file
> wdat1 <- wdat0 %>%
+   filter(!is.na(pair)) %>%
+   group_by(pair) %>%
+   mutate(perc_more = rank(vm.community.subj) - 1) %>%
+   ungroup() %>%
+   mutate(pairF = factor(pair)) %>%
+   select(
+     perc_more, da_prop_vm_20pct_06, vm_change, pair, pairF, social.capital01, community.resp01,
+     vm.community.norm2, vm.community.subj, vcid, dauid, csd_pop_06, csd_prop_vm_20pct_06,
+     da_pop_dens_06, csd_urban_06, csd_pop_dens_06, community_area_km
+   ) %>%
+   droplevels()
> 
> xb1_ranked0 <- xBalance(perc_more ~ da_prop_vm_20pct_06 + vm_change,
+   strata = list(unstr = NULL, pair = ~pair), report = "all", data = wdat1
+ )
> xb1_ranked0$overall
      chisquare df p.value
unstr 0.0001235  2  0.9999
pair  0.3623576  2  0.8343
> xb1_ranked0$results
, , strata = unstr

                     stat
vars                   Control Treatment   adj.diff adj.diff.null.sd   std.diff          z      p
  da_prop_vm_20pct_06  0.07691   0.07691 -1.415e-06         0.003781 -1.218e-05 -0.0003742 0.9997
  vm_change           -0.03037  -0.03034  2.455e-05         0.002225  3.592e-04  0.0110314 0.9912

, , strata = pair

                     stat
vars                   Control Treatment   adj.diff adj.diff.null.sd   std.diff        z      p
  da_prop_vm_20pct_06  0.07691   0.07691 -1.415e-06        2.364e-06 -1.218e-05 -0.59854 0.5495
  vm_change           -0.03037  -0.03034  2.455e-05        3.206e-04  3.592e-04  0.07657 0.9390

attr(,"originals")
[1] "da_prop_vm_20pct_06" "vm_change"          
> 
> xb1_ranked <- balanceTest(perc_more ~ da_prop_vm_20pct_06 + vm_change + strata(pair), data = wdat1, p.adjust = "none")
> xb1_ranked$overall
     chisquare df p.value
pair 0.3624149  2  0.8343
--   0.0001235  2  0.9999
> xb1_ranked$results[, c("Control", "Treatment", "adj.diff", "std.diff", "p"), ] # , "pair"]
, , strata = pair

                     stat
vars                   Control Treatment   adj.diff   std.diff      p
  da_prop_vm_20pct_06  0.07691   0.07691 -1.415e-06 -1.218e-05 0.5494
  vm_change           -0.03037  -0.03034  2.455e-05  3.592e-04 0.9390

, , strata = --

                     stat
vars                   Control Treatment   adj.diff   std.diff      p
  da_prop_vm_20pct_06  0.07691   0.07691 -1.415e-06 -1.218e-05 0.9997
  vm_change           -0.03037  -0.03034  2.455e-05  3.592e-04 0.9912

> 
> 
> ## An alternative for the overall test which should be fairly close to it
> library(formula.tools)
> library(coin)
> coin_fmla <- da_prop_vm_20pct_06 + vm_change ~  vm.community.subj | pairF
> coin_test <- independence_test(coin_fmla, data = wdat1, teststat = "quadratic")
> coin_test_perm <- independence_test(coin_fmla, data = wdat1, teststat = "quadratic", distribution = approximate(nresample = 1000))
> coin::pvalue(coin_test)
[1] 0.4375
> coin::pvalue(coin_test_perm)
[1] 0.463
99 percent confidence interval:
 0.4221 0.5042 

> 
> coin_fmla_rank <- da_prop_vm_20pct_06 + vm_change ~ perc_more | pairF
> coin_test_rank <- independence_test(coin_fmla_rank, data = wdat1, teststat = "quadratic")
> coin_test_rank_perm <- independence_test(coin_fmla_rank,
+   data = wdat1,
+   teststat = "quadratic", distribution = approximate(nresample = 1000)
+ )
> coin::pvalue(coin_test_rank)
[1] 0.8343
> coin::pvalue(coin_test_rank_perm)
[1] 0.826
99 percent confidence interval:
 0.7931 0.8558 

> 
> pairdiffs <- wdat1 %>%
+   group_by(pair) %>%
+   summarize(
+     da_prop_vm_20pct_06 = abs(diff(da_prop_vm_20pct_06)),
+     vm_change = abs(diff(vm_change)),
+     vm.community.subj = abs(diff(vm.community.subj)),
+     csd_pop_06 = abs(diff(csd_pop_06)),
+     csd_pop_dens_06 = abs(diff(csd_pop_dens_06)),
+     csd_prop_vm_20pct_06 = abs(diff(csd_prop_vm_20pct_06)),
+     # da_pop_dens_06 = abs(diff(da_pop_dens_06)),
+     community.area.km = abs(diff(community_area_km))
+   )
> 
> sapply(pairdiffs, summary)
          pair da_prop_vm_20pct_06 vm_change vm.community.subj csd_pop_06 csd_pop_dens_06 csd_prop_vm_20pct_06 community.area.km
Min.       1.0           0.000e+00  0.000000            0.0100        0.0           0.000            0.0000000         1.675e-03
1st Qu.  472.2           0.000e+00  0.000000            0.0700      861.5           2.343            0.0007737         9.951e+00
Median   943.5           0.000e+00  0.007808            0.1700    31284.0          84.224            0.0179847         7.464e+01
Mean     943.5           4.632e-05  0.009957            0.2426   288585.8         459.556            0.0591439         8.224e+02
3rd Qu. 1414.8           4.411e-05  0.018504            0.3300   320419.8         351.662            0.0663463         3.879e+02
Max.    1886.0           3.962e-04  0.030000            4.1900  2501910.0        4552.788            0.5422422         3.662e+04
> 
> sapply(pairdiffs, function(x) {
+   quantile(x, sort(unique(c(seq(0, 1, .1), seq(.9, 1, .01)))), na.rm = TRUE)
+ })
       pair da_prop_vm_20pct_06 vm_change vm.community.subj csd_pop_06 csd_pop_dens_06 csd_prop_vm_20pct_06 community.area.km
0%      1.0           0.0000000  0.000000            0.0100          0            0.00             0.000000         1.675e-03
10%   189.5           0.0000000  0.000000            0.0300          0            0.00             0.000000         1.937e+00
20%   378.0           0.0000000  0.000000            0.0600          0            0.00             0.000000         5.850e+00
30%   566.5           0.0000000  0.000000            0.0800       3476           10.55             0.003316         1.572e+01
40%   755.0           0.0000000  0.003436            0.1200      10457           34.02             0.009075         3.540e+01
50%   943.5           0.0000000  0.007808            0.1700      31284           84.22             0.017985         7.464e+01
60%  1132.0           0.0000000  0.012185            0.2200      72679          160.56             0.030590         1.394e+02
70%  1320.5           0.0000000  0.016475            0.2900     150722          270.55             0.052144         2.665e+02
80%  1509.0           0.0000993  0.020390            0.3800     433140          505.74             0.095595         5.994e+02
90%  1697.5           0.0001970  0.024477            0.5700     811542         1642.19             0.191771         1.607e+03
91%  1716.4           0.0002062  0.025003            0.5835     890718         2338.33             0.198537         1.801e+03
92%  1735.2           0.0002204  0.025628            0.6100     983953         2612.74             0.224051         2.045e+03
93%  1754.1           0.0002324  0.026221            0.6400    1512561         2622.17             0.234841         2.367e+03
94%  1772.9           0.0002555  0.026669            0.6700    1542636         2756.71             0.259034         2.876e+03
95%  1791.8           0.0002701  0.027273            0.7000    1834732         2904.36             0.304956         3.749e+03
96%  1810.6           0.0002913  0.027909            0.7800    1925240         3085.83             0.319108         4.549e+03
97%  1829.5           0.0003117  0.028365            0.8500    2298613         3350.71             0.333054         6.182e+03
98%  1848.3           0.0003376  0.029047            0.9300    2395483         3528.41             0.373991         9.292e+03
99%  1867.2           0.0003612  0.029502            1.1115    2445965         3779.94             0.417343         1.496e+04
100% 1886.0           0.0003962  0.030000            4.1900    2501910         4552.79             0.542242         3.662e+04
> 
> ### facts for the body.tex
> 
> # number of people in final design
> nrow(wdat1)
[1] 3772
> # number of DAs
> length(unique(wdat1$dauid))
[1] 3098
> 
> 
> ### The plot
> 
> ## A plot of the success of the matchings to balance da_prop_vm_20pct_06 but allow variation in vm.community.norm2
> pdf(file = here("Figures_Tables", "nbmplot_anyDA_new.pdf"), width = 6, height = 6)
> par(oma = rep(0, 4), mgp = c(1.5, .5, 0), mar = c(3, 3, 2, 0))
> nbmplot(wdat1,
+   yvar = "da_prop_vm_20pct_06", xvar = "vm.community.norm2", strata = "pair", points = FALSE,
+   xlab = "Perceptions of Local Community Drawings", ylab = "Census Figures",
+   main = "Proportion Visible Minority"
+ )
> abline(0, 1, col = gray(.5), lwd = 1)
> dev.off()
null device 
          1 
> 
> system("touch Data/match_assess_anyDA_new.done")
> 
